﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using MathNet.Numerics;

namespace LPLib
{
    public class TreeAnalyst
    {
        private NucleotideEvolutionModel model;

        private Tree<Polymer<double[]>> thisTree;

        private int length;
        private NucleotidePMF encoding = new NucleotidePMF();
        private GapProbabilities gapProbabilities;
        public GapProbabilities GapProbabilities
        {
            get { return gapProbabilities; }
            set { gapProbabilities = value; }
        }

        private double logMinTime = -7;

        public TreeAnalyst(NucleotideEvolutionModel _model /*, double[] _mutationRates = null */)
        {
            model = _model;
            gapProbabilities = new GapProbabilities();
        }

        private static int sequenceLength(Tree<Polymer<double[]>> tree)
        {
            int length = -1;
            if (tree.Children.Count == 0)
            {
                if (tree.Contents == null)
                {
                    throw new ArgumentException("The tree has an empty leaf.");
                }
                return tree.Contents.Seq.Length;
            }
            else
            {
                foreach (Tree<Polymer<double[]>> child in tree.Children)
                {
                    try
                    {
                        length = sequenceLength(child, length);
                    }
                    catch (ArgumentException e)
                    {
                        throw e;
                    }
                }
                return length;
            }
        }

        private static int sequenceLength(Tree<Polymer<double[]>> tree, int length)
        {
            if (tree.Children.Count == 0)
            {
                if (tree.Contents == null)
                {
                    throw new ArgumentException("The tree has an empty leaf.");
                }

                if (length > -1 && tree.Contents.Seq.Length != length)
                {
                    throw new ArgumentException("Sequences are not all the same length.");
                }

                return tree.Contents.Seq.Length;
            }
            else
            {
                foreach (Tree<Polymer<double[]>> child in tree.Children)
                {
                    try
                    {
                        length = sequenceLength(child, length);
                    }
                    catch (ArgumentException e)
                    {
                        throw e;
                    }
                }
                return length;
            }
        }

        public void ComputePosteriors(Tree<Polymer<double[]>> tree, Polymer<double[]> uca, int len = -1)
        {
            if (len < 0)
            {
                length = sequenceLength(tree);
            }

            // trim or extend uca if its length differs from those of the tree tips
            if (uca.Seq.Length < length)
            {
                // extend with gaps
                int ucaSeqLength = uca.Seq.Length;
                double[][] temp = new double[length][];
                Array.Copy(uca.Seq, temp, uca.Seq.Length);
                uca.Seq = temp;
                for (int i = ucaSeqLength; i < temp.Length; i++)
                {
                    uca.Seq[i] = new double[5];
                    uca.Seq[i][4] = 1.0;
                }            
            }
            else if (uca.Seq.Length > length)
            {
                // trim
                uca.Seq = uca.Subsequence(0, length - 1);
            }

            tree.Contents = uca;
            computePosteriors(tree);
        }

        /// <summary>
        /// Computes the likelihood at every position in the alignment associated with the tree.
        /// </summary>
        /// <param name="tree">The tree being examined.</param>
        /// <param name="likelihood">The array of likelihoods for each nucleotide at each position.</param>
        /// <param name="cdr3S">The position where CDR3 starts.</param>
        public void ComputeLikelihood(Tree<Polymer<double[]>> tree, int len = -1)
        {
            if (len < 0)
            {
                length = sequenceLength(tree);
                //mutationRates = new double[length];
                //TBKMath.Utilities.FillArray(ref mutationRates, 1);
            }

            foreach (Tree<Polymer<double[]>> child in tree.Children)
            {
                ComputeLikelihood(child, length);
            }
            if (tree.Children.Count == 0)
                return;

            double[][] pmf = new double[length][];

            for (int iPos = 0; iPos < length; iPos++)
            {
                pmf[iPos] = new double[NucleotidePMF.Dim];
                for (int iChar = 0; iChar < 4; iChar++)
                {
                    pmf[iPos][iChar] = 1;
                    foreach (Tree<Polymer<double[]>> child in tree.Children)
                    {
                        double p = 0;
                        for (int iCharChild = 0; iCharChild < 4; iCharChild++)
                        {
                            p += child.Contents.Seq[iPos][iCharChild] * model.TransitionProbability(iChar, iCharChild, child.BranchLength);
                        }
                        p += gapProbabilities.OpenDeletionProbability * (1 - model.TransitionProbability(iChar, iChar, child.BranchLength)) * child.Contents.Seq[iPos][4];
                        pmf[iPos][iChar] *= p;
                    }
                }

                pmf[iPos][4] = 1;    
                foreach (Tree<Polymer<double[]>> child in tree.Children)
                {
                    double sum = 0;
                    for (int iChar = 0; iChar < 4; iChar++ )
                    {
                        sum += child.Contents.Seq[iPos][iChar];
                    }
                    // we take the "raw" mutation frequency to be 1 - P(0|0)
                    pmf[iPos][4] *= (gapProbabilities.OpenInsertionProbability * (1 - model.TransitionProbability(0, 0, child.BranchLength)) * sum + child.Contents.Seq[iPos][4] );
                }
            } 
            Polymer<double[]> contents = new Polymer<double[]>(encoding, tree.Name, pmf);
            tree.Contents = contents;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="tree"></param>
        /// <param name="likelihood"></param>
        /// <param name="posteriors"></param>
        private void computePosteriors(Tree<Polymer<double[]>> tree)
        {
            NucleotidePMF pmfEncoding = new NucleotidePMF();

            if (tree.Children.Count == 0) // The tree is a tip. The sequence in Contents is already the posterior.
                return;

            if (tree.Parent != null)
            {
                if (tree.BranchLength == 0)
                {
                    // if the branch length is zero, the posterior for this node is the same as that for the parent
                    double[][] seq = new double[tree.Parent.Contents.Seq.Length][];
                    for (int i = 0; i < seq.Length; i++)
                    {
                        seq[i] = (double[])tree.Parent.Contents.Seq[i].Clone();
                    }
                    Polymer<double[]> pSeq = new Polymer<double[]>(pmfEncoding, tree.Contents.Name, seq);
                    tree.Contents = pSeq;
                }
                else
                {
                    int length = sequenceLength(tree);
                    Polymer<double[]> pSeq = new Polymer<double[]>(pmfEncoding, tree.Contents.Name);
                    pSeq.Seq = new double[length][];
                    for (int iPos = 0; iPos < length; iPos++)
                    {
                        // compute kernel
                        double norm = 0;
                        double[,] kernel = new double[NucleotidePMF.Dim, NucleotidePMF.Dim];
                        for (int parentNuc = 0; parentNuc < 4; parentNuc++)
                        {
                            norm = 0;
                            for (int childNuc = 0; childNuc < 4; childNuc++)
                            {
                                norm += tree.Contents.Seq[iPos][childNuc] * model.TransitionProbability(parentNuc, childNuc, tree.BranchLength);
                            }
                            norm += gapProbabilities.OpenDeletionProbability * (1 - model.TransitionProbability(parentNuc, parentNuc, tree.BranchLength)) * tree.Contents.Seq[iPos][4];
                            for (int childNuc = 0; childNuc < 4; childNuc++)
                            {
                                kernel[childNuc, parentNuc] = tree.Contents.Seq[iPos][childNuc] * model.TransitionProbability(parentNuc, childNuc, tree.BranchLength) / norm;
                            }
                            kernel[4, parentNuc] = gapProbabilities.OpenDeletionProbability * (1 - model.TransitionProbability(parentNuc, parentNuc, tree.BranchLength)) * tree.Contents.Seq[iPos][4] / norm;
                        }

                        // kernel[childNuc,4] 
                        norm = 0;
                        double sum = 0;
                        for (int childNuc = 0; childNuc < 4; childNuc++)
                        {
                            sum += tree.Contents.Seq[iPos][childNuc];
                        }
                        norm += (gapProbabilities.OpenInsertionProbability * (1 - model.TransitionProbability(0, 0, tree.BranchLength)) * sum + tree.Contents.Seq[iPos][4]);
                        for (int childNuc = 0; childNuc < 4; childNuc++)
                        {
                            kernel[childNuc, 4] = gapProbabilities.OpenInsertionProbability * (1 - model.TransitionProbability(0, 0, tree.BranchLength)) * tree.Contents.Seq[iPos][childNuc] / norm;
                        }
                        kernel[4, 4] = tree.Contents.Seq[iPos][4] / norm;

                        double[] p = new double[NucleotidePMF.Dim];
                        for (int childNuc = 0; childNuc < NucleotidePMF.Dim; childNuc++)
                        {
                            for (int parentNuc = 0; parentNuc < NucleotidePMF.Dim; parentNuc++)
                            {
                                p[childNuc] += kernel[childNuc, parentNuc] * tree.Parent.Contents.Seq[iPos][parentNuc];
                            }
                        }
                        pSeq.Seq[iPos] = p;
                    }
                    tree.Contents = pSeq;
                }
            }

            foreach (Tree<Polymer<double[]>> child in tree.Children)
            {
                computePosteriors(child);
            }
        }

        public static void CollectMutations(Tree<Polymer<byte>> tree, Dictionary<string, MutationList> mutations)
        {
            if (tree.Parent != null && tree.Contents != null)
            {
                PolymerPair<byte, byte> pair = new PolymerPair<byte, byte>(tree.Contents, tree.Parent.Contents);
                MutationList mlist = MutationFinder.GetMutations(pair);
            }

            foreach (Tree<Polymer<byte>> t in tree.Children)
            {
                CollectMutations(t, mutations);
            }

        }

        public static double[] TotalFluxes(Tree<Polymer<double[]>> tree, int len = -1)
        {
            if (tree.Parent == null)
                return null;

            if (len < 0)
            {
                len = sequenceLength(tree);
            }

            // loop over sites
            double[] fluxes = new double[4];
            for (int i = 0; i < len; i++)
            {
                // get the fluxes for the nucleotide pair
                fluxes[0] +=  NucleotidePMF.IDFlux(tree.Parent.Contents.Seq[i], tree.Contents.Seq[i]);
                fluxes[1] +=  NucleotidePMF.TransitionFlux(tree.Parent.Contents.Seq[i], tree.Contents.Seq[i]);
                fluxes[2] += NucleotidePMF.TransversionFlux(tree.Parent.Contents.Seq[i], tree.Contents.Seq[i]);
                fluxes[3] += NucleotidePMF.GapFlux(tree.Parent.Contents.Seq[i], tree.Contents.Seq[i]);
            }
            return fluxes;
        }

        private double dLogLDt(double logT)
        {
            // calculates the derivative of the log likelihood with respect to the evolution time 
            // for the branch of the tree leading back to the parent
            if (thisTree.Parent == null)
                return double.NaN;

            double t = Math.Exp(logT);

            // loop over sites
            double dlogldt = 0;
            for (int i = 0; i < length; i++)
            {
                // get the fluxes for the nucleotide pair
                double phiID = NucleotidePMF.IDFlux(thisTree.Parent.Contents.Seq[i], thisTree.Contents.Seq[i]);
                double phiS = NucleotidePMF.TransitionFlux(thisTree.Parent.Contents.Seq[i], thisTree.Contents.Seq[i]);
                double phiV = 1 - phiID - phiS;

                // compute the probabilities for the relevant state changes
                double[] typeProbs = model.TypeProbabilities(t /* * mutationRates[i] */);
                double delta = (phiS / typeProbs[1] - phiID / typeProbs[0]) * (model.Parameters[0]*(1-2*typeProbs[1] - typeProbs[2]) - 2* typeProbs[1] + typeProbs[2])/2;
                delta += (phiV/typeProbs[2]  - phiID/typeProbs[0]) * (1- 2*typeProbs[2]);
                double x = model.Parameters[0];
                // delta *= mutationRates[i];
                dlogldt += delta;
            }

            return dlogldt;
        }

        public double LogL(Tree<Polymer<double[]>> tree, double t)
        {
            // calculates the log likelihood along the branch of the tree leading back to the parent
            if (tree.Parent == null)
                return double.NaN;

            length = tree.Contents.Seq.Length;
            //mutationRates = new double[length];
            //TBKMath.Utilities.FillArray(ref mutationRates, 1);

            // loop over sites
            double logl = 0;
            for (int i = 0; i < length; i++)
            {
                // get the fluxes for the nucleotide pair
                double phiID = NucleotidePMF.IDFlux(tree.Parent.Contents.Seq[i], tree.Contents.Seq[i]);
                double phiS = NucleotidePMF.TransitionFlux(tree.Parent.Contents.Seq[i], tree.Contents.Seq[i]);
                double phiV = 1 - phiID - phiS;

                // compute the probabilities for the relevant state changes
                double[] typeProbs = model.TypeProbabilities(t /* * mutationRates[i] */);

                double delta = phiID * Math.Log(typeProbs[0]) + phiS * Math.Log(typeProbs[1]) + phiV * Math.Log(typeProbs[2]);
                logl += delta;
            }
            return logl;
        }

        public double MLTime(Tree<Polymer<double[]>> tree, double tol = 1.0E-4)
        {
            thisTree = tree;
            length = tree.Contents.Seq.Length;
            double[] fluxes = TotalFluxes(tree);
            
            if (fluxes[0] >= length) // it may be larger by floating-point arithmetic error
                return Math.Exp(logMinTime);

            double muFreq = 1 - fluxes[0] / length;
            
            double tHat = -0.75 * Math.Log(1 - 4.0 / 3 * muFreq);
            double lower = Math.Log(tHat) - 0.69;
            double upper = lower + 1.39;

            TBKMath.Bisector bis = new TBKMath.Bisector(dLogLDt);
            bis.func = new TBKMath.Bisector.FunctionDelegate(dLogLDt);
            double root = bis.getRoot(lower, upper);
            return Math.Exp(root);
        }

        public void ComputeAllMLTimes(Tree<Polymer<double[]>> tree, double tol = 1.0E-4)
        {
            if (!Tree<Polymer<double[]>>.AllNodesOccupied(tree))
                throw new ArgumentException("There are unoccupied nodes in the tree.");

            if (tree.Parent != null)
            {
                tree.BranchLength = MLTime(tree, tol);
            }
            foreach (Tree<Polymer<double[]>> child in tree.Children)
            {
                ComputeAllMLTimes(child);
            }
        }

        public static void TransferGaps(int[][] index, Tree<Polymer<double[]>> tree)
        {
            foreach (Tree<Polymer<double[]>> child in tree.Children)
            {
                TransferGaps(index, child);
            }
            if (tree.Parent != null)
            {
                double[][] newSeq = new double[index.Length][];
                for (int k = 0; k < index.Length; k++)
                {
                    if (index[k][0] > -1)
                    {
                        newSeq[k] = tree.Contents.Seq[index[k][0]];
                    }
                    else
                    {
                        newSeq[k] = new double[] { 0, 0, 0, 0, 1 };
                    }
                }
                tree.Contents.Seq = newSeq;
            }
        }
    }

    public class TreeAnalystWithIndex
    {
        private NucleotideEvolutionModel model;
        private GapScores gapScores;
        private Tree<Polymer<double[]>> thisTree;
        private NucleotidePMF encoding = new NucleotidePMF();
        public AlignmentIndex aIndex;
        private double logDistMin = -7;

        public TreeAnalystWithIndex(NucleotideEvolutionModel _model, GapScores gapScores, AlignmentIndex _index)
        {
            model = _model;
            aIndex = _index;
            this.gapScores = gapScores;
        }

        public void ComputePosteriors(Tree<Polymer<double[]>> tree, Polymer<double[]> uca, int len = -1)
        {
            // trim or extend uca if its length differs from those of the tree tips
            if (uca.Seq.Length < aIndex.Length)
            {
                // extend with gaps
                int ucaSeqLength = uca.Seq.Length;
                double[][] temp = new double[aIndex.Length][];
                Array.Copy(uca.Seq, temp, uca.Seq.Length);
                uca.Seq = temp;
                for (int i = ucaSeqLength; i < temp.Length; i++)
                {
                    uca.Seq[i] = new double[5];
                    uca.Seq[i][4] = 1.0;
                }
            }
            else if (uca.Seq.Length > aIndex.Length)
            {
                // trim
                uca.Seq = uca.Subsequence(0, aIndex.Length - 1);
            }

            tree.Contents = uca;
            computePosteriors(tree);
        }

        /// <summary>
        /// Computes the likelihood at every position in the alignment associated with the tree.
        /// </summary>
        /// <param name="tree">The tree being examined.</param>
        /// <param name="likelihood">The array of likelihoods for each nucleotide at each position.</param>
        /// <param name="cdr3S">The position where CDR3 starts.</param>
        public void ComputeLikelihoodSequence(Tree<Polymer<double[]>> tree)
        {
            foreach (Tree<Polymer<double[]>> child in tree.Children)
            {
                ComputeLikelihoodSequence(child);
            }
            if (tree.Children.Count == 0)
                return;

            double[][] pmf = new double[aIndex.Length][];
            int[] newIndex = new int[aIndex.Length];

            for (int iPos = 0; iPos < aIndex.Length; iPos++)
            {
                pmf[iPos] = new double[NucleotidePMF.Dim];
                for (int iChar = 0; iChar < NucleotidePMF.Dim; iChar++)
                {
                    pmf[iPos][iChar] = 1;
                    foreach (Tree<Polymer<double[]>> child in tree.Children)
                    {
                        int i = aIndex[child.Contents.Name,iPos];
                        if (i < 0)
                        {
                            continue;
                        }

                        double p = 0;
                        for (int iCharChild = 0; iCharChild < NucleotidePMF.Dim; iCharChild++)
                        {
                            p += child.Contents.Seq[i][iCharChild] * model.TransitionProbability(iChar, iCharChild, child.BranchLength /* * mutationRates[iPos] */);
                        }
                        pmf[iPos][iChar] *= p;
                    }
                }
            }

            Polymer<double[]> contents = new Polymer<double[]>(encoding, tree.Name, pmf);
            tree.Contents = contents;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="tree"></param>
        /// <param name="likelihood"></param>
        /// <param name="posteriors"></param>
        private void computePosteriors(Tree<Polymer<double[]>> tree)
        {
            NucleotidePMF pmfEncoding = new NucleotidePMF();

            if (tree.Children.Count == 0) // The tree is a tip. The sequence in Contents is already the posterior.
                return;

            if (tree.Parent != null)
            {
                if (tree.BranchLength == 0)
                {
                    // if the branch length is zero, the posterior for this node is the same as that for the parent
                    double[][] seq = new double[tree.Parent.Contents.Seq.Length][];
                    for (int i = 0; i < seq.Length; i++)
                    {
                        seq[i] = (double[])tree.Parent.Contents.Seq[i].Clone();
                    }
                    Polymer<double[]> pSeq = new Polymer<double[]>(pmfEncoding, tree.Contents.Name, seq);
                    tree.Contents = pSeq;
                }
                else
                {
                    Polymer<double[]> pSeq = new Polymer<double[]>(pmfEncoding, tree.Contents.Name);
                    pSeq.Seq = new double[aIndex.Length][];
                    for (int iPos = 0; iPos < aIndex.Length; iPos++)
                    {
                        // compute kernel
                        double[,] kernel = new double[NucleotidePMF.Dim, NucleotidePMF.Dim];
                        for (int parentNuc = 0; parentNuc < NucleotidePMF.Dim; parentNuc++)
                        {
                            double norm = 0;
                            for (int nuc = 0; nuc < NucleotidePMF.Dim; nuc++)
                            {
                                norm += tree.Contents.Seq[iPos][nuc] * model.TransitionProbability(parentNuc, nuc, tree.BranchLength /* * mutationRates[iPos] */);
                            }
                            for (int iChar = 0; iChar < NucleotidePMF.Dim; iChar++)
                            {
                                kernel[iChar, parentNuc] = tree.Contents.Seq[iPos][iChar] * model.TransitionProbability(parentNuc, iChar, tree.BranchLength /* * mutationRates[iPos] */) / norm;
                            }
                        }

                        double[] p = new double[NucleotidePMF.Dim];
                        for (int nuc = 0; nuc < NucleotidePMF.Dim; nuc++)
                        {
                            for (int parentNuc = 0; parentNuc < NucleotidePMF.Dim; parentNuc++)
                            {
                                p[nuc] += kernel[nuc, parentNuc] * tree.Parent.Contents.Seq[iPos][parentNuc];
                            }
                        }
                        pSeq.Seq[iPos] = p;
                    }
                    tree.Contents = pSeq;
                }
            }

            foreach (Tree<Polymer<double[]>> child in tree.Children)
            {
                computePosteriors(child);
            }
        }

        public static void CollectMutations(Tree<Polymer<byte>> tree, Dictionary<string, MutationList> mutations)
        {
            if (tree.Parent != null && tree.Contents != null)
            {
                PolymerPair<byte, byte> pair = new PolymerPair<byte, byte>(tree.Contents, tree.Parent.Contents);
                MutationList mlist = MutationFinder.GetMutations(pair);
            }

            foreach (Tree<Polymer<byte>> t in tree.Children)
            {
                CollectMutations(t, mutations);
            }

        }

        public static double[] TotalFluxes(Tree<Polymer<double[]>> tree, AlignmentIndex aindex)
        {
            if (tree.Parent == null)
                return null;

            // loop over sites
            double[] fluxes = new double[5];
            for (int i = 0; i < aindex.Length; i++)
            {
                int iPosParent = aindex[tree.Parent.Contents.Name, i];
                int iPosSelf = aindex[tree.Contents.Name, i];
                if ( iPosSelf > -1 && iPosParent > -1)
                {
                    // get the fluxes for the nucleotide pair
                    fluxes[0] += NucleotidePMF.IDFlux(tree.Parent.Contents.Seq[iPosParent], tree.Contents.Seq[iPosSelf]);
                    fluxes[1] += NucleotidePMF.TransitionFlux(tree.Parent.Contents.Seq[iPosParent], tree.Contents.Seq[iPosSelf]);
                    fluxes[2] += NucleotidePMF.TransversionFlux(tree.Parent.Contents.Seq[iPosParent], tree.Contents.Seq[iPosSelf]);
                }
                else if (iPosSelf < 0) // deletion
                {
                    fluxes[3]++;
                }
                else // insertion
                {
                    fluxes[4]++;
                }
            }
            return fluxes;
        }

        private double dLogLDt(double logT)
        {
            // calculates the derivative of the log likelihood with respect to the evolution time 
            // for the branch of the tree leading back to the parent
            if (thisTree.Parent == null)
                return double.NaN;

            double t = Math.Exp(logT);

            // loop over sites
            double dlogldt = 0;
            for (int i = 0; i < aIndex.Length; i++)
            {
                double phiID = 0;
                double phiS = 0;
                int iParent = aIndex[thisTree.Parent.Contents.Name, i];
                int iChild = aIndex[thisTree.Contents.Name, i];
                if (iParent > -1 &&  iChild> -1)
                {
                    // get the fluxes for the nucleotide pair
                    phiID = NucleotidePMF.IDFlux(thisTree.Parent.Contents.Seq[iParent], thisTree.Contents.Seq[iChild]);
                    phiS = NucleotidePMF.TransitionFlux(thisTree.Parent.Contents.Seq[iParent],
                        thisTree.Contents.Seq[iChild]);
                }
                double phiV = 1 - phiID - phiS;

                // compute the probabilities for the relevant state changes
                double[] typeProbs = model.TypeProbabilities(t /* * mutationRates[i] */);
                double delta = (phiS / typeProbs[1] - phiID / typeProbs[0]) * (model.Parameters[0] * (1 - 2 * typeProbs[1] - typeProbs[2]) - 2 * typeProbs[1] + typeProbs[2]) / 2;
                delta += (phiV / typeProbs[2] - phiID / typeProbs[0]) * (1 - 2 * typeProbs[2]);
                double x = model.Parameters[0];
                // delta *= mutationRates[i];
                dlogldt += delta;
            }

            return dlogldt;
        }

        public double LogL(Tree<Polymer<double[]>> tree, double t)
        {
            // calculates the log likelihood along the branch of the tree leading back to the parent
            if (tree.Parent == null)
                return double.NaN;

            // loop over sites
            double logl = 0;
            for (int i = 0; i < aIndex.Length; i++)
            {
                double phiID = 0;
                double phiS = 0;
                int iParent = aIndex[thisTree.Parent.Contents.Name, i];
                int iChild = aIndex[thisTree.Contents.Name, i];
                if (iParent > -1 && iChild > -1)
                {
                    // get the fluxes for the nucleotide pair
                    phiID = NucleotidePMF.IDFlux(thisTree.Parent.Contents.Seq[iParent], thisTree.Contents.Seq[iChild]);
                    phiS = NucleotidePMF.TransitionFlux(thisTree.Parent.Contents.Seq[iParent],
                        thisTree.Contents.Seq[iChild]);
                }
                double phiV = 1 - phiID - phiS;

                // compute the probabilities for the relevant state changes
                double[] typeProbs = model.TypeProbabilities(t /* * mutationRates[i] */);

                double delta = phiID * Math.Log(typeProbs[0]) + phiS * Math.Log(typeProbs[1]) + phiV * Math.Log(typeProbs[2]);
                logl += delta;
            }
            return logl;
        }

        public double MLTime(Tree<Polymer<double[]>> tree, double tol = 1.0E-4)
        {
            thisTree = tree;
            double[] fluxes = TotalFluxes(tree, aIndex);

            if (fluxes[0] >= aIndex.Length) // there are no mutations. (Use GTE because it may be larger by floating-point arithmetic error.)
                return Math.Exp(logDistMin);

            double muFreq = 1 - fluxes[0] / aIndex.Length;

            double tHat = -0.75 * Math.Log(1 - 4.0 / 3 * muFreq);
            double lower = Math.Log(tHat) - 0.69;
            double upper = lower + 1.39;

            TBKMath.Bisector bis = new TBKMath.Bisector(dLogLDt);
            bis.func = new TBKMath.Bisector.FunctionDelegate(dLogLDt);
            double root = bis.getRoot(lower, upper);
            return Math.Exp(root);
        }

        public void ComputeAllMLTimes(Tree<Polymer<double[]>> tree, double tol = 1.0E-4)
        {
            if (!Tree<Polymer<double[]>>.AllNodesOccupied(tree))
                throw new ArgumentException("There are unoccupied nodes in the tree.");

            if (tree.Parent != null)
            {
                tree.BranchLength = MLTime(tree, tol);
            }
            foreach (Tree<Polymer<double[]>> child in tree.Children)
            {
                ComputeAllMLTimes(child);
            }
        }
    }

    public class TrefoilAnalyst
    {
        private NucleotideEvolutionModel model;
        private GapProbabilities gapProbabilities;
        private Trefoil trefoil;
        private NucleotidePMF encoding = new NucleotidePMF();
        private double minTime = 0.001;
        private int activeMember;
        double tolerance = 1.0E-4;

        public int ActiveMember
        {
            get { return activeMember; }
            set { activeMember = value; }
        }

        public TrefoilAnalyst(NucleotideEvolutionModel model, GapProbabilities gapProbabilities, Trefoil trefoil)
        {
            this.model = model;
            this.trefoil = trefoil;
            this.gapProbabilities = gapProbabilities;
        }

        /// <summary>
        /// Computes the likelihood at every position in the alignment associated with the tree.
        /// </summary>
        /// <param name="tree">The tree being examined.</param>
        /// <param name="likelihood">The array of likelihoods for each nucleotide at each position.</param>
        /// <param name="cdr3S">The position where CDR3 starts.</param>
        public void ComputeLikelihoodSequencesSegmented()
        {
            List<double[]> rootSequence = new List<double[]>();
            List<double[]> intermediateSequence = new List<double[]>();
            double indelProbabilities = 1;
            for (int segment = 0; segment < trefoil.NSegments; segment++ )
            {
                if (trefoil.SegmentStart[segment][3] < 0) // this segment is not present in the root sequence
                {
                    double pIndel;
                    if ((trefoil.SegmentStart[segment][0] > -1 && trefoil.SegmentStart[segment][1] > -1)) /* both of the terminal polymers contain the segment
                        so we must account for descent from the MRCA */
                    {

                        pIndel = 1 - Math.Exp(-gapProbabilities.OpenInsertionProbability* trefoil.DistanceToRoot);

                        // the segment in the MRCA begins at its current length
                        trefoil.SegmentStart[segment][2] = intermediateSequence.Count;
                        for (int position = 0; position < trefoil.SegmentLength[segment]; position++)
                        {
                            double[] pmf = new double[] { 1, 1, 1, 1, 0 };
                            backPropagate(pmf, trefoil.Polymers[0].Seq[trefoil.SegmentStart[segment][0] + position], trefoil.DistancesToMRCA[0]);
                            backPropagate(pmf, trefoil.Polymers[1].Seq[trefoil.SegmentStart[segment][1] + position], trefoil.DistancesToMRCA[1]);
                            intermediateSequence.Add(pmf);
                            
                            double pThisPosition = 0;
                            for (int nuc = 0; nuc < 4; nuc++)
                            {
                                pThisPosition += pmf[nuc];
                            }
                            indelProbabilities *= pThisPosition;
                        }
                    }
                    else if (trefoil.SegmentStart[segment][1] > -1)
                    {
                        // the insertion occurs in the branch from the MRCA to polymer 1
                        pIndel = 1 - Math.Exp(-gapProbabilities.OpenInsertionProbability * trefoil.DistancesToMRCA[1]);
                    }
                    else // note that the segment cannot be absent in both polymers
                    {
                        // the insertion occurs in the branch from the MRCA to polymer 0
                        pIndel = pIndel = 1 - Math.Exp(gapProbabilities.OpenInsertionProbability * trefoil.DistancesToMRCA[1]);
                    }

                    indelProbabilities *= Math.Pow(gapProbabilities.ContinueInsertionProbability/4, trefoil.SegmentLength[segment] - 1);
                }
                else // the segment is present in the root
                {
                    double mu;
                    if (trefoil.SegmentStart[segment][0] > -1 && trefoil.SegmentStart[segment][1] > -1) /* both of the terminal polymers contains the segment
                        so we must account for descent from the MRCA */
                    {
                        // the segment in the MRCA begins at its current length
                        trefoil.SegmentStart[segment][2] = intermediateSequence.Count;
                        for (int position = 0; position < trefoil.SegmentLength[segment]; position++)
                        {
                            double[] ancestorPMF = new double[] { 1, 1, 1, 1, 0 };
                            backPropagate(ancestorPMF, trefoil.Polymers[0].Seq[trefoil.SegmentStart[segment][0] + position], trefoil.DistancesToMRCA[0]);
                            backPropagate(ancestorPMF, trefoil.Polymers[1].Seq[trefoil.SegmentStart[segment][1] + position], trefoil.DistancesToMRCA[1]);
                            double[] rootPMF = new double[] {1, 1, 1 ,1 , 0};
                            intermediateSequence.Add(ancestorPMF);
                            backPropagate(rootPMF, ancestorPMF, trefoil.DistanceToRoot);
                            rootSequence.Add(rootPMF);
                        }
                    }
                    else if (trefoil.SegmentStart[segment][0] > -1 || trefoil.SegmentStart[segment][1] > -1)
                    {
                        // the segment in the MRCA begins at its current length
                        trefoil.SegmentStart[segment][2] = intermediateSequence.Count;
                        // determine which polymer experienced the deletion
                        if (trefoil.SegmentStart[segment][0] < 0)
                        {
                            // the deletion occurs in the branch from the MRCA to polymer 0
                            mu = 1 - Math.Exp(-gapProbabilities.OpenDeletionProbability * trefoil.DistancesToMRCA[0]);
                        }
                        else
                        {
                            // the deletion occurs in the branch from the MRCA to polymer 1
                            mu = 1 - Math.Exp(-gapProbabilities.OpenDeletionProbability * trefoil.DistancesToMRCA[1]);
                        }

                        indelProbabilities *= mu * Math.Pow(gapProbabilities.ContinueDeletionProbability, trefoil.SegmentLength[segment] - 1);
                        
                        for (int position = 0; position < trefoil.SegmentLength[segment]; position++)
                        {
                            double[] ancestorPMF = new double[] { 1, 1, 1, 1, 0 };
                            if (trefoil.SegmentStart[segment][0] > -1)
                                backPropagate(ancestorPMF, trefoil.Polymers[0].Seq[trefoil.SegmentStart[segment][0] + position], trefoil.DistancesToMRCA[0]);

                            if (trefoil.SegmentStart[segment][1] > -1)
                                backPropagate(ancestorPMF, trefoil.Polymers[1].Seq[trefoil.SegmentStart[segment][1] + position], trefoil.DistancesToMRCA[1]);

                            double[] rootPMF = new double[] { 1, 1, 1, 1, 0 };
                            intermediateSequence.Add(ancestorPMF);
                            backPropagate(rootPMF, ancestorPMF, trefoil.DistanceToRoot);
                            rootSequence.Add(rootPMF);
                        }
                    }
                    else
                    {
                        // the segment exists only in the root
                        if (segment!=0 && segment != trefoil.NSegments - 1) // don't charge for leading or trailing root nucleotides
                        {
                            // the deletion occurs in the branch from the root to the MRCA
                            for (int position = 0; position < trefoil.SegmentLength[segment]; position++)
                            {
                                rootSequence.Add(new double[] { 0.25, 0.25, 0.25, 0.25, 0 });
                            }
                            mu = 1 - model.TransitionProbability(1, 1, trefoil.DistanceToRoot);
                            indelProbabilities += mu * Math.Pow(gapProbabilities.ContinueDeletionProbability, trefoil.SegmentLength[segment] - 1);
                        }
                    }
                }
            }

            trefoil.MRCALikelihoodSequence = new Polymer<double[]>(encoding, "Intermediate", intermediateSequence.ToArray());
            trefoil.RootLikelihoodSequence = new Polymer<double[]>(encoding, "UCA", rootSequence.ToArray());
            trefoil.IndelPartialLikelihood = indelProbabilities;
        }

        private void backPropagate(double[] pmf, double[] child, double time)
        {
            for (int iNuc = 0; iNuc < 4; iNuc++)
            {
                double p = 0;
                for (int cNuc = 0; cNuc < 4; cNuc++)
                {
                    p += child[cNuc] * model.TransitionProbability(iNuc, cNuc, time);
                }
                pmf[iNuc] *= p;
            }
        }

        private void dBackPropagateDT(double[] pmf, double[] child, double time)
        {
            for (int iNuc = 0; iNuc < 4; iNuc++)
            {
                double p = 0;
                for (int cNuc = 0; cNuc < 4; cNuc++)
                {
                    p += child[cNuc] * model.dTransitionProbabilityDT(iNuc, cNuc, time);
                }
                pmf[iNuc] *= p;
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="tree"></param>
        /// <param name="likelihood"></param>
        /// <param name="posteriors"></param>
        //private void computePosteriors()
        //{
        //    // prior on root
        //    if (trefoil.Distances[2] == 0)
        //    {
        //        // if the branch length is zero, the posterior for this node is the same as that for the parent
        //        Polymer<double[]> pSeq = new Polymer<double[]>(encoding, trefoil.IntermediateSequence.Name, trefoil.RootSequence.Seq); 
        //    }
        //    else
        //    {
        //        Polymer<double[]> pSeq = new Polymer<double[]>(encoding, trefoil.IntermediateSequence.Name);
        //        List<double[]> seq = new List<double[]>();
        //        for (int iPos = 0; iPos < trefoil.AlignmentLength; iPos++)
        //        {
        //            // compute kernel
        //            double[,] kernel = new double[NucleotidePMF.Dim, NucleotidePMF.Dim];
        //            for (int parentNuc = 0; parentNuc < NucleotidePMF.Dim; parentNuc++)
        //            {
        //                double norm = 0;
        //                for (int nuc = 0; nuc < NucleotidePMF.Dim; nuc++)
        //                {
        //                    norm += trefoil.IntermediateSequence.Seq[iPos][nuc] * model.TransitionProbability(parentNuc, nuc, tree.BranchLength);
        //                }
        //                for (int iChar = 0; iChar < NucleotidePMF.Dim; iChar++)
        //                {
        //                    kernel[iChar, parentNuc] = tree.Contents.Seq[iPos][iChar] * model.TransitionProbability(parentNuc, iChar, tree.BranchLength /* * mutationRates[iPos] */) / norm;
        //                }
        //            }

        //            double[] p = new double[NucleotidePMF.Dim];
        //            for (int nuc = 0; nuc < NucleotidePMF.Dim; nuc++)
        //            {
        //                for (int parentNuc = 0; parentNuc < NucleotidePMF.Dim; parentNuc++)
        //                {
        //                    p[nuc] += kernel[nuc, parentNuc] * tree.Parent.Contents.Seq[iPos][parentNuc];
        //                }
        //            }
        //            pSeq.Seq[iPos] = p;
        //        }
        //        tree.Contents = pSeq;
        //    }
        //}

        public double[] TotalFluxes()
        {
            double[] fluxes = new double[4];
            // loop over sites
            if (activeMember < 2) // fluxes from an observed sequence to the intermediate
            {
                for (int k = 0; k < trefoil.AlignmentLength; k++)
                {
                    int iPosParent = trefoil.Index[k][2];
                    int iPosSelf = trefoil.Index[k][activeMember];
                    if (iPosSelf > -1 && iPosParent > -1)
                    {
                        // get the fluxes for the nucleotide pair
                        fluxes[0] += NucleotidePMF.IDFlux(trefoil.MRCANormalized.Seq[iPosParent], trefoil.Polymers[activeMember].Seq[iPosSelf]);
                        fluxes[1] += NucleotidePMF.TransitionFlux(trefoil.MRCANormalized.Seq[iPosParent], trefoil.Polymers[activeMember].Seq[iPosSelf]);
                        fluxes[2] += NucleotidePMF.TransversionFlux(trefoil.MRCANormalized.Seq[iPosParent], trefoil.Polymers[activeMember].Seq[iPosSelf]);
                    }
                    else if (iPosSelf > -1 || iPosParent > -1)
                    {
                        fluxes[3] += 1;
                    }
                }
            }
            else
            {
                for (int k = 0; k < trefoil.AlignmentLength; k++)
                {
                    int iPosParent = trefoil.Index[k][3];
                    int iPosSelf = trefoil.Index[k][activeMember];
                    if (iPosSelf > -1 && iPosParent > -1)
                    {
                        // get the fluxes for the nucleotide pair
                        fluxes[0] += NucleotidePMF.IDFlux(trefoil.RootNormalized.Seq[iPosParent], trefoil.MRCANormalized.Seq[iPosSelf]);
                        fluxes[1] += NucleotidePMF.TransitionFlux(trefoil.RootNormalized.Seq[iPosParent], trefoil.MRCANormalized.Seq[iPosSelf]);
                        fluxes[2] += NucleotidePMF.TransversionFlux(trefoil.RootNormalized.Seq[iPosParent], trefoil.MRCANormalized.Seq[iPosSelf]);
                    }
                    else if (iPosSelf > -1 || iPosParent > -1)
                    {
                        fluxes[3] += 1;
                    }
                }
            }
            return fluxes;
        }

        public double[] DirectFluxesBetweenObserved()
        {
            double[] fluxes = new double[4];
            // loop over sites
            for (int k = 0; k < trefoil.AlignmentLength; k++)
            {
                int iPos1 = trefoil.Index[k][0];
                int iPos2 = trefoil.Index[k][1];
                if (iPos1 > -1 && iPos2 > -1)
                {
                    // get the fluxes for the nucleotide pair
                    fluxes[0] += NucleotidePMF.IDFlux(trefoil.Polymers[0].Seq[iPos1], trefoil.Polymers[1].Seq[iPos2]);
                    fluxes[1] += NucleotidePMF.TransitionFlux(trefoil.Polymers[0].Seq[iPos1], trefoil.Polymers[1].Seq[iPos2]);
                    fluxes[2] += NucleotidePMF.TransversionFlux(trefoil.Polymers[0].Seq[iPos1], trefoil.Polymers[1].Seq[iPos2]);
                }
                else
                {
                    // gap flux
                    fluxes[3] += 1;
                }
            }
            return fluxes;
        }

        private double dLogLDt(double t)
        {
            // calculates the derivative of the log likelihood with respect to the evolution time 

            if (activeMember > 2)
            {
                return double.NaN;
            }

            double dLogLDt = 0;
            for (int segment = 0; segment < trefoil.NSegments; segment++)
            {
                // check for indels
                if (trefoil.SegmentStart[segment][3] > -1) // segment present in root 
                {
                    if (trefoil.SegmentStart[segment][2] > -1) // segment present in MRCA
                    {                                           // must be present in at least one polymer
                        if (trefoil.SegmentStart[segment][1] > -1) // segment present in Polymer0
                        {
                            if (trefoil.SegmentStart[segment][0] < 0) // segment not present in all segments
                            {
                                if (activeMember == 0)
                                // deletion along branch 0
                                {
                                    dLogLDt += gapProbabilities.OpenDeletionProbability
                                    / (Math.Exp(gapProbabilities.OpenDeletionProbability * t) - 1);
                                }
                            }
                            // no indel if present in all polymers
                        }
                        else
                        {
                            // deletion along branch 1
                            if (activeMember == 1)
                            {
                                dLogLDt += gapProbabilities.OpenDeletionProbability
                                        / (Math.Exp(gapProbabilities.OpenDeletionProbability * t) - 1);
                            }
                        }
                    }
                    else // segment not present in MCRA
                    {
                        if (activeMember == 2)
                        {
                            dLogLDt += gapProbabilities.OpenDeletionProbability
                                / (Math.Exp(gapProbabilities.OpenDeletionProbability * t) - 1);
                        }
                    }
                }
                else
                {
                    if (trefoil.SegmentStart[segment][2] > -1) // segment present in MRCA
                    {                                           // must be present in both polymers
                        // insertion in MRCA
                        if (activeMember == 2)
                        {
                            dLogLDt += gapProbabilities.OpenInsertionProbability
                                / (Math.Exp(gapProbabilities.OpenInsertionProbability * t) - 1);
                        }
                    }
                    else // segment not present in MCRA, must be present in one polymer
                    {
                        if (trefoil.SegmentStart[segment][1] > -1) // segment present in Polymer0
                        {
                            // insertion in polymer 0
                            if (activeMember == 0)
                            {
                                dLogLDt += gapProbabilities.OpenInsertionProbability
                                    / (Math.Exp(gapProbabilities.OpenInsertionProbability * t) - 1);
                            }
                        }
                        else
                        {
                            // insertion in polymer 1
                            if (activeMember == 1)
                            {
                                dLogLDt += gapProbabilities.OpenInsertionProbability
                                    / (Math.Exp(gapProbabilities.OpenInsertionProbability * t) - 1);
                            }
                        }
                    }
                }

                double[] lA;
                double[] dLA;
                double[] lR;
                double[] dLR;

                for (int i = 0; i < trefoil.SegmentLength[segment]; i++)
                {
                    lA = new double[] { 1, 1, 1, 1, 0 };
                    dLA = new double[] { 1, 1, 1, 1, 0 };

                    if (trefoil.SegmentStart[segment][0] > -1)
                    {
                        if (activeMember == 0)
                        {
                            backPropagate(lA, trefoil.Polymers[0].Seq[trefoil.SegmentStart[segment][0] + i], t);
                            dBackPropagateDT(dLA, trefoil.Polymers[0].Seq[trefoil.SegmentStart[segment][0] + i], t);
                        }
                        else
                        {
                            backPropagate(lA, trefoil.Polymers[0].Seq[trefoil.SegmentStart[segment][0] + i], trefoil.DistancesToMRCA[0]);
                        }
                    }

                    if (trefoil.SegmentStart[segment][1] > -1)
                    {
                        if (activeMember == 1)
                        {
                            backPropagate(lA, trefoil.Polymers[1].Seq[trefoil.SegmentStart[segment][1] + i], t);
                            dBackPropagateDT(dLA, trefoil.Polymers[1].Seq[trefoil.SegmentStart[segment][1] + i], t);
                        }
                        else
                        {
                            backPropagate(lA, trefoil.Polymers[1].Seq[trefoil.SegmentStart[segment][1] + i], trefoil.DistancesToMRCA[1]);
                        }
                    }

                    // for derivative
                    if (trefoil.SegmentStart[segment][0] > -1 && activeMember == 1)
                    {
                        backPropagate(dLA, trefoil.Polymers[0].Seq[trefoil.SegmentStart[segment][0] + i], trefoil.DistancesToMRCA[0]);
                    }
                    else if (trefoil.SegmentStart[segment][1] > -1 && activeMember == 0)
                    {
                        backPropagate(dLA, trefoil.Polymers[1].Seq[trefoil.SegmentStart[segment][1] + i], trefoil.DistancesToMRCA[1]);
                    }

                    // for likelihood
                    if (trefoil.SegmentStart[segment][3] > -1)
                    {
                        lR = new double[] { 1, 1, 1, 1, 0 };
                        dLR = new double[] { 1, 1, 1, 1, 0 };
                        if (activeMember != 2)
                        {
                            backPropagate(lR, lA, trefoil.DistanceToRoot);
                            backPropagate(dLR, dLA, trefoil.DistanceToRoot);
                        }
                        else
                        {
                            backPropagate(lR, lA, t);
                            dBackPropagateDT(dLR, lA, t);
                        }
                        double numerator = 0;
                        double denominator = 0;
                        for (int r = 0; r < 4; r++)
                        {
                            numerator += trefoil.RootNormalized.Seq[trefoil.SegmentStart[segment][3] + i][r] * dLR[r];
                            denominator += trefoil.RootNormalized.Seq[trefoil.SegmentStart[segment][3] + i][r] * lR[r];
                        }
                        dLogLDt += numerator / denominator;
                    }
                }
            }

            return dLogLDt;
        }

        private double[] dLogLDts(double[] u)
        {
            // calculates the derivative of the log likelihood with respect to the evolution time 
            // for the branch of the tree leading back to the parent

            double[] t = new double[3];
            for (int i = 0; i < 3; i++)
            {
                t[i] = Math.Max(u[i],0);
            }

            double[] dLogLIndel = new double[3];
            double[] dLogLDt = new double[3];
            for (int segment = 0; segment < trefoil.NSegments; segment++)
            {
                bool[] insertion = new bool[] { false, false, false };
                // check for indels
                if (trefoil.SegmentStart[segment][3] > -1) // segment present in root 
                {
                    if (trefoil.SegmentStart[segment][2] > -1) // segment present in MRCA
                    {                                           // must be present in at least one polymer
                        if (trefoil.SegmentStart[segment][1] > -1) // segment present in Polymer0
                        {
                            if (trefoil.SegmentStart[segment][0] < 0) // segment not present in all segments
                            {
                                // deletion along branch 0
                                dLogLIndel[0] = gapProbabilities.OpenDeletionProbability 
                                    / (Math.Exp(gapProbabilities.OpenDeletionProbability * t[0]) - 1)
                                    + (trefoil.SegmentLength[segment] - 1) * Math.Log(gapProbabilities.ContinueDeletionProbability);
                            }
                            // no indel if present in all polymers
                        }
                        else
                        {
                            // deletion along branch 1
                            dLogLIndel[1] = gapProbabilities.OpenDeletionProbability
                                / (Math.Exp(gapProbabilities.OpenDeletionProbability * t[1]) - 1)
                                + (trefoil.SegmentLength[segment] - 1) * Math.Log(gapProbabilities.ContinueDeletionProbability);
                        }
                    }
                    else // segment not present in MCRA
                    {
                        dLogLIndel[2] = gapProbabilities.OpenDeletionProbability
                            / (Math.Exp(gapProbabilities.OpenDeletionProbability * t[2]) - 1)
                            + (trefoil.SegmentLength[segment] - 1) * Math.Log(gapProbabilities.ContinueDeletionProbability);
                    }
                }
                else
                {
                    if (trefoil.SegmentStart[segment][2] > -1) // segment present in MRCA
                    {                                           // must be present in both polymers
                        // insertion in MRCA
                        dLogLIndel[2] = gapProbabilities.OpenInsertionProbability
                            / (Math.Exp(gapProbabilities.OpenInsertionProbability * t[2]) - 1);
                    }
                    else // segment not present in MCRA, must be present in one polymer
                    {
                        if (trefoil.SegmentStart[segment][1] > -1) // segment present in Polymer0
                        {
                            // insertion in polymer 0
                            dLogLIndel[0] = gapProbabilities.OpenInsertionProbability
                                / (Math.Exp(gapProbabilities.OpenInsertionProbability * t[0]) - 1);
                        }
                        else
                        {
                            // insertion in polymer 1
                            dLogLIndel[1] = gapProbabilities.OpenInsertionProbability
                                / (Math.Exp(gapProbabilities.OpenInsertionProbability * t[1]) - 1);
                        }
                    }
                }

                double[] lA;
                double[] dLA0;
                double[] dLA1;
                double[] lR;
                double[] dLR0;
                double[] dLR1;
                double[] dLRt;

                for (int i = 0; i < trefoil.SegmentLength[segment]; i++)
                {
                    lA = new double[] { 1, 1, 1, 1, 0 };
                    dLA0 = new double[] { 1, 1, 1, 1, 0 };
                    dLA1 = new double[] { 1, 1, 1, 1, 0 };

                    if (trefoil.SegmentStart[segment][0] > -1)
                    {
                        backPropagate(lA, trefoil.Polymers[0].Seq[trefoil.SegmentStart[segment][0] + i], t[0]);
                        dBackPropagateDT(dLA0, trefoil.Polymers[0].Seq[trefoil.SegmentStart[segment][0] + i], t[0]);
                    }
                    else
                    {
                    }

                    if (trefoil.SegmentStart[segment][1] > -1)
                    {
                        backPropagate(lA, trefoil.Polymers[1].Seq[trefoil.SegmentStart[segment][1] + i], t[1]);
                        dBackPropagateDT(dLA1, trefoil.Polymers[1].Seq[trefoil.SegmentStart[segment][1] + i], t[1]);
                    }
                    else
                    {
                        // get indel component of derivative
                    }

                    // for derivative
                    if (trefoil.SegmentStart[segment][0] > -1)
                    {
                        backPropagate(dLA1, trefoil.Polymers[0].Seq[trefoil.SegmentStart[segment][0] + i], t[0]);
                    }
                    else
                    {

                    }

                    if (trefoil.SegmentStart[segment][1] > -1)
                    {
                        backPropagate(dLA0, trefoil.Polymers[1].Seq[trefoil.SegmentStart[segment][1] + i], t[1]);
                    }
                    else
                    {

                    }

                    // for likelihood
                    if (trefoil.SegmentStart[segment][3] > -1)
                    {
                        lR = new double[] { 1, 1, 1, 1, 0 };
                        backPropagate(lR, lA, t[2]);
                        dLR0 = new double[] { 1, 1, 1, 1, 0 };
                        dLR1 = new double[] { 1, 1, 1, 1, 0 };
                        dLRt = new double[] { 1, 1, 1, 1, 0 };
                        backPropagate(dLR0, dLA0, t[2]);
                        backPropagate(dLR1, dLA1, t[2]);
                        dBackPropagateDT(dLRt, lA, t[2]);
                        double numerator0 = 0;
                        double numerator1 = 0;
                        double numeratort = 0;
                        double denominator = 0;
                        for (int r = 0; r < 4; r++)
                        {
                            numerator0 += trefoil.RootNormalized.Seq[trefoil.SegmentStart[segment][3] + i][r] * dLR0[r];
                            numerator1 += trefoil.RootNormalized.Seq[trefoil.SegmentStart[segment][3] + i][r] * dLR1[r];
                            numeratort += trefoil.RootNormalized.Seq[trefoil.SegmentStart[segment][3] + i][r] * dLRt[r];
                            denominator += trefoil.RootNormalized.Seq[trefoil.SegmentStart[segment][3] + i][r] * lR[r];
                        }
                        dLogLDt[0] += numerator0 / denominator;
                        dLogLDt[1] += numerator1 / denominator;
                        dLogLDt[2] += numeratort / denominator;
                    }
                    else
                    {
                        
                    }
                }
            }

            double alpha = Math.Pow(trefoil.Index.Length, 3);
            for (int i = 0; i < 3; i++)
            {
                if (u[i] < 0 && dLogLDt[i] < 0)
                {
                    dLogLDt[i] = u[i] * u[i];
                }
            }
            return dLogLDt;
        }

        public double LogL(double t)
        {
            // calculates the log likelihood
            // for the branch of the tree leading back to the parent

            // loop over sites
            double logl = 0;
            for (int k = 0; k < trefoil.AlignmentLength; k++)
            {
                double phiID = 0;
                double phiS = 0;
                double phiV = 1;
                int iChild = trefoil.Index[k][activeMember];
                int iParent;
                if (activeMember < 2)
                {
                    iParent = trefoil.Index[k][activeMember];
                    if (iParent > -1 && iChild > -1)
                    {
                        // get the fluxes for the nucleotide pair
                        phiID = NucleotidePMF.IDFlux(trefoil.MRCANormalized.Seq[iParent], trefoil.Polymers[activeMember].Seq[iChild]);
                        phiS = NucleotidePMF.TransitionFlux(trefoil.MRCANormalized.Seq[iParent], trefoil.Polymers[activeMember].Seq[iChild]);
                    }
                }
                else if (activeMember == 2)
                 {
                    iParent = trefoil.Index[k][activeMember];
                    if (iParent > -1 && iChild > -1)
                    {
                        // get the fluxes for the nucleotide pair
                        phiID = NucleotidePMF.IDFlux(trefoil.RootNormalized.Seq[iParent], trefoil.MRCANormalized.Seq[iChild]);
                        phiS = NucleotidePMF.TransitionFlux(trefoil.RootNormalized.Seq[iParent], trefoil.MRCANormalized.Seq[iChild]);
                    }
                }
                phiV = 1 - phiID - phiS;

                // compute the probabilities for the relevant state changes
                double[] typeProbs = model.TypeProbabilities(t);
                double delta = phiID * Math.Log(typeProbs[0]) + phiS * Math.Log(typeProbs[1]) + phiV * Math.Log(typeProbs[2]);
                logl += delta;
            }
            return logl;
        }

        public void ComputeMLTimes()
        {
            // initialize the times to a reasonable value. There may be more reasonable values to use, but we'll see.
            double tReasonable = 0.05;
            double[] t = new double[]{tReasonable, tReasonable, tReasonable};
            trefoil.DistancesToMRCA[0] = t[0];
            trefoil.DistancesToMRCA[1] = t[1];
            trefoil.DistanceToRoot = t[2];
            double[] tLast = new double[3];
            double[] lower = new double[]{0,0,0};
            double[] upper = new double[3]{tReasonable, tReasonable, tReasonable};
            double error = double.MaxValue;
            // initialize the upper bounds
            for (activeMember = 0; activeMember < 3; activeMember++ )
            {
                if (dLogLDt(minTime) < 0)
                {
                    continue;
                }
                else
                {
                    while (dLogLDt(upper[activeMember]) > 0)
                    {
                        upper[activeMember] *= 1.5;
                    }
                }
            }
            while (error > tolerance)
            {
                error = 0;
                for (activeMember = 0; activeMember < 3; activeMember++)
                {
                    tLast[activeMember] = t[activeMember];
                    if (dLogLDt(minTime) < 0)
                    {
                        t[activeMember] = minTime;
                        lower[activeMember] = minTime;
                        upper[activeMember] = 0.01;
                    }
                    else
                    {
                        while (dLogLDt(lower[activeMember]) < 0)
                        {
                            lower[activeMember] *= 0.9;
                        }
                        while (dLogLDt(upper[activeMember]) > 0)
                        {
                            upper[activeMember] *= 1.1;
                        }
                        t[activeMember] = MathNet.Numerics.FindRoots.OfFunction(dLogLDt, lower[activeMember], upper[activeMember]);
                        lower[activeMember] = t[activeMember];
                        upper[activeMember] = t[activeMember];
                    }
                    error += Math.Pow(t[activeMember] - tLast[activeMember], 2);
                }
                trefoil.DistancesToMRCA[0] = t[0];
                trefoil.DistancesToMRCA[1] = t[1];
                trefoil.DistanceToRoot = t[2];
            }
        }

    }

    public static class PairAnalyst
    {
        public static double[] TotalFluxes(PolymerPair<double[], double[]> pair)
        {
            double[] fluxes = new double[4];
            // loop over sites
            for (int k = 0; k < pair.Index.Length; k++)
            {
                int iPosParent = pair.Index[k][1];
                int iPosSelf = pair.Index[k][0];
                if (iPosSelf > -1 && iPosParent > -1)
                {
                    // get the fluxes for the nucleotide pair
                    fluxes[0] += NucleotidePMF.IDFlux(pair.Polymer1.Seq[iPosParent], pair.Polymer0.Seq[iPosSelf]);
                    fluxes[1] += NucleotidePMF.TransitionFlux(pair.Polymer1.Seq[iPosParent], pair.Polymer0.Seq[iPosSelf]);
                    fluxes[2] += NucleotidePMF.TransversionFlux(pair.Polymer1.Seq[iPosParent], pair.Polymer0.Seq[iPosSelf]);
                }
                else
                {
                    fluxes[3] += 1;
                }
            }
            return fluxes;
        }

        public static Polymer<double[]> GetConsensus(PolymerPair<double[], double[]> pair)
        {
            Polymer<double[]> consensus = new Polymer<double[]>(new NucleotidePMF(), "consensus");
            consensus.Seq = new double[pair.Index.Length][];
            for (int i = 0; i < pair.Index.Length; i++)
            {
                consensus.Seq[i] = new double[]{1,1,1,1,1};
                double sum = 0;
                for (int nuc = 0; nuc < 5; nuc++)
                {
                    if (pair.Index[i][0] > -1)
                        consensus.Seq[i][nuc] *= pair.Polymer0.Seq[pair.Index[i][0]][nuc];

                    if (pair.Index[i][1] > -1)
                        consensus.Seq[i][nuc] *= pair.Polymer1.Seq[pair.Index[i][1]][nuc];

                    sum += consensus.Seq[i][nuc];
                }
                if (sum > 0)
                {
                    for (int nuc = 0; nuc < 5; nuc++)
                    {
                        consensus.Seq[i][nuc] /= sum;
                    }
                }
                else
                {
                    for (int nuc = 0; nuc < 5; nuc++)
                    {
                        consensus.Seq[i][nuc] = 0.5*(pair.Polymer0.Seq[pair.Index[i][0]][nuc] + pair.Polymer1.Seq[pair.Index[i][1]][nuc]);
                    }
                }
            }
            return consensus;
        }
    }

    public class AlignmentIndex
    {
        private int length = 0;
        public int Length
        {
            get { return length; }
        }

        private Dictionary<string, int[]> index;

        public AlignmentIndex(string name, int[] index)
        {
            this.index = new Dictionary<string, int[]>();
            this.index.Add(name, index);
            length = index.Length;
        }

        public int this[string name, int position]
        {
            get { return index[name][position]; }
            set { index[name][position] = value;}
        }

        public void AddMember(string name)
        {
            index.Add(name, new int[length]);
        }

        public void AddMember(string name, int[] newIndex)
        {
            if (newIndex.Length != length && length > 0)
            {
                throw new ArgumentException("The new index must be the same length as existing indices.");
            }
            if (length == 0)
            {
                length = newIndex.Length;
            }

            index.Add(name, newIndex);
        }

        public void RemoveMember(string name)
        {
            index.Remove(name);
        }
    }
}
